
NUM_MC = 10000
set.seed(0)

# SIMULATION PARAMETERS ################################################
N = 75
T = 75
b = matrix(c(0.15,0.2,0.15),,1)
DGP_AR = TRUE ### AR if TRUE; MA if FALSE ###
rho = .25 ### In case of AR ###
q = 10
MAcoef = matrix(c(rep(1/q,q)),,1) ### In case of MA ###
true_beta = matrix(c(.1,.1),2,1)

# LIST OF ALTERNATIVE SHIFTS ###########################################
list_alt = (-50:50)/100

# TUNING FOR THOMPSON ##################################################
m_thompson = 2

# MENZEL BOOTSTRAP ITERATIONS ##########################################
num_boot = 1000
# MNW BOOTSTRAP ITERATIONS #############################################
num_boot_MNW = 1000

# MAMEN RANDOM VARIABLE ################################################
mammen <- function(para){ a = para[1]; b = para[2]; p = para[3]; return( (p*a+(1-p)*b)^2 + (p*a^2+(1-p)*b^2-1)^2 + (p*a^3+(1-p)*b^3-1)^2 ); }
mammen_para = nlm(mammen,c(-1,1,0.5))$estimate
rmammen <- function(n,para){(para[1]-para[2])*(runif(n)<para[3])+para[2]}

# FUNCTION: EIGENVALUE-CORRECTION ######################################
EV_Correction <- function( A ){
 EV = eigen(A)
 Eval = EV$values
 Eval = (Eval >= 0) * Eval + (Eval < 0 ) * 0
 Evec = EV$vectors
 return( (Evec) %*% diag(Eval) %*% t(Evec) )
}

########################################################################
# BEGIN MONTE CARLO ####################################################
########################################################################
list_IID_95 = NULL
list_PNW_95 = NULL
list_CRi_95 = NULL
list_CRt_95 = NULL
list_CGM_95 = NULL
list_MNW_95 = NULL
list_MEN_95 = NULL
list_ME2_95 = NULL
list_THO_95 = NULL
list_DDG_95 = NULL
list_CHS_95 = NULL
for( iter in 1:NUM_MC ){

########################################################################
# DATA GENERATION (BEGIN)
########################################################################

### ALPHA (NORMAL) (N by 1) ############################################
alpha_X = matrix(rnorm(N,0,1),N,1)
alpha_U = matrix(rnorm(N,0,1),N,1)

### GAMMA (AR NORMAL) (1 by T) #########################################
if( DGP_AR ){
 gamma_X = matrix(rnorm(1,0,1),1,T)
 gamma_U = matrix(rnorm(1,0,1),1,T)
 for( idx in 2:T ){
  gamma_X[1,idx] = rho*gamma_X[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
  gamma_U[1,idx] = rho*gamma_U[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
 }
}

### GAMMA (MA NORMAL) (1 by T) #########################################
if( !DGP_AR ){
 pre_gamma_X = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 pre_gamma_U = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 gamma_X = matrix(0,1,T)
 gamma_U = matrix(0,1,T)
 for( idx in 1:T ){
  gamma_X[1,idx] = t(MAcoef) %*% matrix(pre_gamma_X[1:(dim(MAcoef)[1])+idx-1],,1)
  gamma_U[1,idx] = t(MAcoef) %*% matrix(pre_gamma_U[1:(dim(MAcoef)[1])+idx-1],,1)
 }
}

### EPSILON (NORMAL) (N by T) ##########################################
epsilon_X = matrix(rnorm(N*T,0,1),N,T)
epsilon_U = matrix(rnorm(N*T,0,1),N,T)

X = b[1,1] * alpha_X %*% matrix(1,1,T) +
    b[2,1] * matrix(1,N,1) %*% gamma_X +
    b[3,1] * epsilon_X

U = b[1,1] * alpha_U %*% matrix(1,1,T) +
    b[2,1] * matrix(1,N,1) %*% gamma_U +
    b[3,1] * epsilon_U

Y = true_beta[1,1] + true_beta[2,1] * X + U

########################################################################
# DATA GENERATION (ENDED)
########################################################################



########################################################################
# ESTIMATION (BEGIN)
########################################################################

### OLS ################################################################
beta_hat = solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% t(cbind(1,matrix(c(X),,1))) %*% matrix(c(Y),,1)

### RESIDUALS ##########################################################
U_hat = Y - beta_hat[1,1] - beta_hat[2,1] * X

########################################################################
# ESTIMATION (ENDED)
########################################################################



########################################################################
# ESTIMATE IID OMEGA (BEGIN)
########################################################################
Omega_hat_iid = matrix(0,2,2)
Omega_hat_iid[1,1] = t(matrix(c(1*U_hat),,1)) %*% matrix(c(1*U_hat),,1)
Omega_hat_iid[2,1] = t(matrix(c(X*U_hat),,1)) %*% matrix(c(1*U_hat),,1)
Omega_hat_iid[1,2] = t(matrix(c(1*U_hat),,1)) %*% matrix(c(X*U_hat),,1)
Omega_hat_iid[2,2] = t(matrix(c(X*U_hat),,1)) %*% matrix(c(X*U_hat),,1)
Omega_hat_iid = 1/(N*T) * Omega_hat_iid
########################################################################
# ESTIMATE IID OMEGA (ENDED)
########################################################################



########################################################################
# MENZEL'S BOOTSTRAP (BEGIN)
########################################################################

Ubar = mean(U_hat)
XUbar = mean(X*U_hat)
a1 = rowMeans(U_hat)
a2 = rowMeans(X*U_hat)
g1 = colMeans(U_hat)
g2 = colMeans(X*U_hat)
w1 = U_hat - matrix(a1,N,T) - t(matrix(g1,T,N))
w2 = X*U_hat - matrix(a2,N,T) - t(matrix(g2,T,N))

S2_a1 = rowSums((U_hat - Ubar)^2)/(N-1)
S2_a2 = rowSums((X*U_hat - XUbar)^2)/(N-1)
S2_g1 = colSums((U_hat - Ubar)^2)/(T-1)
S2_g2 = colSums((X*U_hat - XUbar)^2)/(T-1)
S2_w1 = sum((U_hat- matrix(a1,N,T) - t(matrix(g1,T,N)) - Ubar)^2) / (N*T - N - T)
S2_w2 = sum((X*U_hat - matrix(a2,N,T) - t(matrix(g2,T,N)) - XUbar)^2) / (N*T - N - T)

sigma2_a1 = max(c(0,S2_a1 - S2_w1/T))
sigma2_a2 = max(c(0,S2_a2 - S2_w2/T))
sigma2_g1 = max(c(0,S2_g1 - S2_w1/N))
sigma2_g2 = max(c(0,S2_g2 - S2_w2/N))
sigma2_w1 = S2_w1
sigma2_w2 = S2_w2

eps_hat = U_hat - matrix(a1,N,T) - t(matrix(g1,T,N)) - mean(U_hat)
mu4_e = sqrt(2)*max(c(0.1,sqrt(mean(eps_hat^4)-mean(eps_hat)^4)))
kappa_a = mu4_e * log(T)
kappa_g = mu4_e * log(N)

D_a1 = 1*(T*sigma2_a1 >= kappa_a)
D_a2 = 1*(T*sigma2_a2 >= kappa_a)
D_g1 = 1*(N*sigma2_g1 >= kappa_g)
D_g2 = 1*(N*sigma2_g2 >= kappa_g)

lambda_a1 = (D_a1 * T * sigma2_a1) / (D_a1*T*sigma2_a1 + sigma2_w1)
lambda_a2 = (D_a2 * T * sigma2_a2) / (D_a2*T*sigma2_a2 + sigma2_w2)
lambda_g1 = (D_g1 * N * sigma2_g1) / (D_g1*N*sigma2_g1 + sigma2_w1)
lambda_g2 = (D_g2 * N * sigma2_g2) / (D_g2*N*sigma2_g2 + sigma2_w2)

lambda_a1_No_Model_Sel = (1 * T * sigma2_a1) / (1*T*sigma2_a1 + sigma2_w1)
lambda_a2_No_Model_Sel = (1 * T * sigma2_a2) / (1*T*sigma2_a2 + sigma2_w2)
lambda_g1_No_Model_Sel = (1 * N * sigma2_g1) / (1*N*sigma2_g1 + sigma2_w1)
lambda_g2_No_Model_Sel = (1 * N * sigma2_g2) / (1*N*sigma2_g2 + sigma2_w2)

boot_beta = NULL
for ( boot_iter in 1:num_boot ){

# DRAW BOOTSTRAP SAMPLES ###############################################
ki = sample(1:N,N,replace=TRUE)
kt = sample(1:T,T,replace=TRUE)

# STEP (a) #############################################################
a1_star = a1[ki]
a2_star = a2[ki]
g1_star = g1[kt]
g2_star = g2[kt]
# STEP (b) #############################################################
omega1_star = matrix(rmammen(N,mammen_para),N,T)
omega2_star = t(matrix(rmammen(T,mammen_para),T,N))
w1_star = omega1_star*omega2_star*w1[ki,kt]
w2_star = omega1_star*omega2_star*w2[ki,kt]
# STEP (c) #############################################################
z1_star = lambda_a1^0.5*matrix(a1_star,N,T) + lambda_g1^0.5*t(matrix(g1_star,T,N)) + w1_star
z2_star = lambda_a2^0.5*matrix(a2_star,N,T) + lambda_g2^0.5*t(matrix(g2_star,T,N)) + w2_star
z1_star_No_Model_Sel = lambda_a1_No_Model_Sel^0.5*matrix(a1_star,N,T) + lambda_g1_No_Model_Sel^0.5*t(matrix(g1_star,T,N)) + w1_star
z2_star_No_Model_Sel = lambda_a2_No_Model_Sel^0.5*matrix(a2_star,N,T) + lambda_g2_No_Model_Sel^0.5*t(matrix(g2_star,T,N)) + w2_star
# STEP (d) #############################################################
boot_beta = cbind(boot_beta, beta_hat + N*T*solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% matrix(c(mean(z1_star),mean(z2_star)),2,1))
boot_beta_No_Model_Sel = cbind(boot_beta, beta_hat + N*T*solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% matrix(c(mean(z1_star_No_Model_Sel),mean(z2_star_No_Model_Sel)),2,1))
}

powers_MEN = NULL
powers_ME2 = NULL

for( idx in 1:length(list_alt) ){
alt = list_alt[idx]
pow_MEN = abs(beta_hat-true_beta-alt) > apply(boot_beta,1,sd)*qnorm(0.975)
pow_ME2 = abs(beta_hat-true_beta-alt) > apply(boot_beta_No_Model_Sel,1,sd)*qnorm(0.975)

powers_MEN = cbind(powers_MEN,pow_MEN[2])
powers_ME2 = cbind(powers_ME2,pow_ME2[2])
}

list_MEN_95 = rbind(list_MEN_95, powers_MEN)
list_ME2_95 = rbind(list_ME2_95, powers_ME2)

########################################################################
# MENZEL'S BOOTSTRAP (ENDED)
########################################################################



########################################################################
# ESTIMATE CHS OMEGA (BEGIN)
########################################################################

### OMEGA HAT 1ST TERM #################################################
Omega_hat_1 = matrix(0,2,2)
for( i in 1:N ){
 Omega_hat_1[1,1] = Omega_hat_1[1,1] + sum( t(matrix(    1*U_hat[i,],1,T)) %*% matrix(    1*U_hat[i,],1,T) )
 Omega_hat_1[2,1] = Omega_hat_1[2,1] + sum( t(matrix(X[i,]*U_hat[i,],1,T)) %*% matrix(    1*U_hat[i,],1,T) )
 Omega_hat_1[1,2] = Omega_hat_1[1,2] + sum( t(matrix(    1*U_hat[i,],1,T)) %*% matrix(X[i,]*U_hat[i,],1,T) )
 Omega_hat_1[2,2] = Omega_hat_1[2,2] + sum( t(matrix(X[i,]*U_hat[i,],1,T)) %*% matrix(X[i,]*U_hat[i,],1,T) )
}
Omega_hat_1 = min(c(N,T))/N/(N*T^2) * Omega_hat_1

### OMEGA HAT 2ND TERM #################################################
Omega_hat_2 = matrix(0,2,2)
for( t in 1:T ){
 Omega_hat_2[1,1] = Omega_hat_2[1,1] + sum( matrix(    1*U_hat[,t],N,1) %*% t(matrix(    1*U_hat[,t],N,1)) )
 Omega_hat_2[2,1] = Omega_hat_2[2,1] + sum( matrix(X[,t]*U_hat[,t],N,1) %*% t(matrix(    1*U_hat[,t],N,1)) )
 Omega_hat_2[1,2] = Omega_hat_2[1,2] + sum( matrix(    1*U_hat[,t],N,1) %*% t(matrix(X[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[2,2] = Omega_hat_2[2,2] + sum( matrix(X[,t]*U_hat[,t],N,1) %*% t(matrix(X[,t]*U_hat[,t],N,1)) )
}
Omega_hat_2 = min(c(N,T))/T/(N^2*T) * Omega_hat_2

### OMEGA HAT 3RD TERM #################################################
Omega_hat_3 = matrix(0,2,2)
Omega_hat_3[1,1] = t(matrix(c(1*U_hat),,1)) %*% matrix(1*c(1*U_hat),,1)
Omega_hat_3[2,1] = t(matrix(c(X*U_hat),,1)) %*% matrix(1*c(1*U_hat),,1)
Omega_hat_3[1,2] = t(matrix(c(1*U_hat),,1)) %*% matrix(1*c(X*U_hat),,1)
Omega_hat_3[2,2] = t(matrix(c(X*U_hat),,1)) %*% matrix(1*c(X*U_hat),,1)
Omega_hat_3 = min(c(N,T))/(N*T)/(N*T) * Omega_hat_3

### OMEGA HAT 4TH TERM FOR THOMPSON ####################################
Omega_hat_4_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_4_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  temp_Omega_hat_4_thompson[1,1] = temp_Omega_hat_4_thompson[1,1] + sum( matrix(      1*U_hat[,t],N,1) %*% t(matrix(      1*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[2,1] = temp_Omega_hat_4_thompson[2,1] + sum( matrix(  X[,t]*U_hat[,t],N,1) %*% t(matrix(      1*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[1,2] = temp_Omega_hat_4_thompson[1,2] + sum( matrix(      1*U_hat[,t],N,1) %*% t(matrix(X[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[2,2] = temp_Omega_hat_4_thompson[2,2] + sum( matrix(  X[,t]*U_hat[,t],N,1) %*% t(matrix(X[,t-j]*U_hat[,t-j],N,1)) )
 }
 temp_Omega_hat_4_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_4_thompson
 Omega_hat_4_thompson = Omega_hat_4_thompson + temp_Omega_hat_4_thompson
}
Omega_hat_4_thompson = min(c(N,T))/T * Omega_hat_4_thompson

### OMEGA HAT 5TH TERM FOR THOMPSON ####################################
Omega_hat_5_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_5_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  temp_Omega_hat_5_thompson[1,1] = temp_Omega_hat_5_thompson[1,1] + sum( matrix(      1*U_hat[,t-j],N,1) %*% t(matrix(      1*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[2,1] = temp_Omega_hat_5_thompson[2,1] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(      1*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[1,2] = temp_Omega_hat_5_thompson[1,2] + sum( matrix(      1*U_hat[,t-j],N,1) %*% t(matrix(  X[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[2,2] = temp_Omega_hat_5_thompson[2,2] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(  X[,t]*U_hat[,t],N,1)) )
 }
 temp_Omega_hat_5_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_5_thompson
 Omega_hat_5_thompson = Omega_hat_5_thompson + temp_Omega_hat_5_thompson
}
Omega_hat_5_thompson = min(c(N,T))/T * Omega_hat_5_thompson

### OMEGA HAT 6TH TERM FOR THOMPSON ####################################
Omega_hat_6_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_6_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  temp_Omega_hat_6_thompson[1,1] = temp_Omega_hat_6_thompson[1,1] + sum( matrix(      1*U_hat[,t],N,1) * (matrix(      1*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[2,1] = temp_Omega_hat_6_thompson[2,1] + sum( matrix(  X[,t]*U_hat[,t],N,1) * (matrix(      1*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[1,2] = temp_Omega_hat_6_thompson[1,2] + sum( matrix(      1*U_hat[,t],N,1) * (matrix(X[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[2,2] = temp_Omega_hat_6_thompson[2,2] + sum( matrix(  X[,t]*U_hat[,t],N,1) * (matrix(X[,t-j]*U_hat[,t-j],N,1)) )
 }
 temp_Omega_hat_6_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_6_thompson
 Omega_hat_6_thompson = Omega_hat_6_thompson + temp_Omega_hat_6_thompson
}
Omega_hat_6_thompson = min(c(N,T))/T * Omega_hat_6_thompson

### OMEGA HAT 5TH TERM FOR THOMPSON ####################################
Omega_hat_7_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_7_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  temp_Omega_hat_7_thompson[1,1] = temp_Omega_hat_7_thompson[1,1] + sum( matrix(      1*U_hat[,t-j],N,1) * (matrix(      1*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[2,1] = temp_Omega_hat_7_thompson[2,1] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) * (matrix(      1*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[1,2] = temp_Omega_hat_7_thompson[1,2] + sum( matrix(      1*U_hat[,t-j],N,1) * (matrix(  X[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[2,2] = temp_Omega_hat_7_thompson[2,2] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) * (matrix(  X[,t]*U_hat[,t],N,1)) )
 }
 temp_Omega_hat_7_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_7_thompson
 Omega_hat_7_thompson = Omega_hat_7_thompson + temp_Omega_hat_7_thompson
}
Omega_hat_7_thompson = min(c(N,T))/T * Omega_hat_7_thompson

### COMPUTE M HAT ######################################################
S = colSums(X*U_hat)[2:T]
Slag = colSums(X*U_hat)[1:(T-1)]
rho_hat = cov(S,Slag)/var(Slag)
m = 1.8171 * (rho_hat^2/(1-rho_hat^2)^2)^(1/3) * T^(1/3)
m = min(c(m,T-1))

### OMEGA HAT 4TH TERM WITH BARTLET KERNEL #############################
Omega_hat_4 = matrix(0,2,2)
if( m >= 1 ){
 for( j in 1:m ){
  temp_Omega_hat_4 = matrix(0,2,2)
  w = 1 - (j/(m+1))
  for( t in (j+1):T ){
   temp_Omega_hat_4[1,1] = temp_Omega_hat_4[1,1] + sum( matrix(      1*U_hat[,t],N,1) %*% t(matrix(      1*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[2,1] = temp_Omega_hat_4[2,1] + sum( matrix(  X[,t]*U_hat[,t],N,1) %*% t(matrix(      1*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[1,2] = temp_Omega_hat_4[1,2] + sum( matrix(      1*U_hat[,t],N,1) %*% t(matrix(X[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[2,2] = temp_Omega_hat_4[2,2] + sum( matrix(  X[,t]*U_hat[,t],N,1) %*% t(matrix(X[,t-j]*U_hat[,t-j],N,1)) )
  }
  temp_Omega_hat_4 = w / (N^2*T) * temp_Omega_hat_4
  Omega_hat_4 = Omega_hat_4 + temp_Omega_hat_4
 }
 Omega_hat_4 = min(c(N,T))/T * Omega_hat_4
}

### OMEGA HAT 5TH TERM WITH BARTLET KERNEL #############################
Omega_hat_5 = matrix(0,2,2)
if( m >= 1 ){
 for( j in 1:m ){
  temp_Omega_hat_5 = matrix(0,2,2)
  w = 1 - (j/(m+1))
  for( t in (j+1):T ){
   temp_Omega_hat_5[1,1] = temp_Omega_hat_5[1,1] + sum( matrix(      1*U_hat[,t-j],N,1) %*% t(matrix(      1*U_hat[,t],N,1)) )
   temp_Omega_hat_5[2,1] = temp_Omega_hat_5[2,1] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(      1*U_hat[,t],N,1)) )
   temp_Omega_hat_5[1,2] = temp_Omega_hat_5[1,2] + sum( matrix(      1*U_hat[,t-j],N,1) %*% t(matrix(  X[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[2,2] = temp_Omega_hat_5[2,2] + sum( matrix(X[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(  X[,t]*U_hat[,t],N,1)) )
  }
  temp_Omega_hat_5 = w / (N^2*T) * temp_Omega_hat_5
  Omega_hat_5 = Omega_hat_5 + temp_Omega_hat_5
 }
 Omega_hat_5 = min(c(N,T))/T * Omega_hat_5
}

########################################################################
# ESTIMATE CHS OMEGA (BEGIN)
########################################################################



########################################################################
# VARIANCE ESTIMATION (BEGIN)
########################################################################

### DESIGN MATRIX Q HAT ################################################
Q_hat = t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) / (N*T)

### STANDARD ERRORS ####################################################
SE_IID = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_iid )                                                                                                                       )%*%solve(Q_hat)/(N*T) ) )
SE_CRi = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 )                                                                                                                         )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CRt = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_2 )                                                                                                                         )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CGM = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 )                                                                                             )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_DDG = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 )                                                                                                           )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_THO = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 + Omega_hat_4_thompson + Omega_hat_5_thompson - Omega_hat_6_thompson - Omega_hat_7_thompson ) )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CHS = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 + Omega_hat_4 + Omega_hat_5 )                                                                 )%*%solve(Q_hat)/min(c(N,T)) ) )

########################################################################
# VARIANCE ESTIMATION (ENDED)
########################################################################



########################################################################
# WILD BOOTSTRAP CRITICAL VALUE (MNW2021) (BEGIN)
########################################################################
MNW_t_list = NULL
for( boot_iterations in 1:num_boot_MNW ){
rademacher = matrix(-1+2*(runif(N*T)<0.5),N,T)
U_hat_wild = rademacher * U_hat
Y_boot = beta_hat[1,1] + beta_hat[2,1] * X + U_hat_wild
beta_hat_boot = solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% t(cbind(1,matrix(c(X),,1))) %*% matrix(c(Y_boot),,1)
U_hat_boot = Y_boot - beta_hat_boot[1,1] - beta_hat_boot[2,1] * X

### OMEGA HAT 1ST TERM #################################################
Omega_hat_1_boot = matrix(0,2,2)
for( i in 1:N ){
 Omega_hat_1_boot[1,1] = Omega_hat_1_boot[1,1] + sum( t(matrix(    1*U_hat_boot[i,],1,T)) %*% matrix(    1*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[2,1] = Omega_hat_1_boot[2,1] + sum( t(matrix(X[i,]*U_hat_boot[i,],1,T)) %*% matrix(    1*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[1,2] = Omega_hat_1_boot[1,2] + sum( t(matrix(    1*U_hat_boot[i,],1,T)) %*% matrix(X[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[2,2] = Omega_hat_1_boot[2,2] + sum( t(matrix(X[i,]*U_hat_boot[i,],1,T)) %*% matrix(X[i,]*U_hat_boot[i,],1,T) )
}
Omega_hat_1_boot = min(c(N,T))/N/(N*T^2) * Omega_hat_1_boot

### OMEGA HAT 2ND TERM #################################################
Omega_hat_2_boot = matrix(0,2,2)
for( t in 1:T ){
 Omega_hat_2_boot[1,1] = Omega_hat_2_boot[1,1] + sum( matrix(    1*U_hat_boot[,t],N,1) %*% t(matrix(    1*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[2,1] = Omega_hat_2_boot[2,1] + sum( matrix(X[,t]*U_hat_boot[,t],N,1) %*% t(matrix(    1*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[1,2] = Omega_hat_2_boot[1,2] + sum( matrix(    1*U_hat_boot[,t],N,1) %*% t(matrix(X[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[2,2] = Omega_hat_2_boot[2,2] + sum( matrix(X[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X[,t]*U_hat_boot[,t],N,1)) )
}
Omega_hat_2_boot = min(c(N,T))/T/(N^2*T) * Omega_hat_2_boot

### OMEGA HAT 3RD TERM #################################################
Omega_hat_3_boot = matrix(0,2,2)
Omega_hat_3_boot[1,1] = t(matrix(c(1*U_hat_boot),,1)) %*% matrix(1*c(1*U_hat_boot),,1)
Omega_hat_3_boot[2,1] = t(matrix(c(X*U_hat_boot),,1)) %*% matrix(1*c(1*U_hat_boot),,1)
Omega_hat_3_boot[1,2] = t(matrix(c(1*U_hat_boot),,1)) %*% matrix(1*c(X*U_hat_boot),,1)
Omega_hat_3_boot[2,2] = t(matrix(c(X*U_hat_boot),,1)) %*% matrix(1*c(X*U_hat_boot),,1)
Omega_hat_3_boot = min(c(N,T))/(N*T)/(N*T) * Omega_hat_3_boot

SE_CGM_boot = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1_boot + Omega_hat_2_boot - Omega_hat_3_boot ) )%*%solve(Q_hat)/min(c(N,T)) ) )
t_wild_boot = (beta_hat_boot - beta_hat ) / matrix(SE_CGM_boot,2,1)
MNW_t_list = cbind( MNW_t_list, t_wild_boot )
}
crit_MNW = apply(abs(MNW_t_list),1,quantile,p=0.95)
########################################################################
# WILD BOOTSTRAP CRITICAL VALUE (MNW2021) (ENDED)
########################################################################

powers_IID_95 = NULL
powers_CRi_95 = NULL
powers_CRt_95 = NULL
powers_CGM_95 = NULL
powers_MNW_95 = NULL
powers_THO_95 = NULL
powers_DDG_95 = NULL
powers_CHS_95 = NULL

for( idx in 1:length(list_alt) ){
alt = list_alt[idx]

pow_IID_95 = t(abs(beta_hat - true_beta - alt) > SE_IID*qnorm(0.975)); if( is.na(pow_IID_95)[1] ){pow_IID_95[1]=0}; if( is.na(pow_IID_95)[2] ){pow_IID_95[2]=0};
pow_CRi_95 = t(abs(beta_hat - true_beta - alt) > SE_CRi*qnorm(0.975)); if( is.na(pow_CRi_95)[1] ){pow_CRi_95[1]=0}; if( is.na(pow_CRi_95)[2] ){pow_CRi_95[2]=0};
pow_CRt_95 = t(abs(beta_hat - true_beta - alt) > SE_CRt*qnorm(0.975)); if( is.na(pow_CRt_95)[1] ){pow_CRt_95[1]=0}; if( is.na(pow_CRt_95)[2] ){pow_CRt_95[2]=0};
pow_CGM_95 = t(abs(beta_hat - true_beta - alt) > SE_CGM*qnorm(0.975)); if( is.na(pow_CGM_95)[1] ){pow_CGM_95[1]=0}; if( is.na(pow_CGM_95)[2] ){pow_CGM_95[2]=0};
pow_MNW_95 = t(abs(beta_hat - true_beta - alt) > SE_CGM*crit_MNW);     if( is.na(pow_CGM_95)[1] ){pow_CGM_95[1]=0}; if( is.na(pow_CGM_95)[2] ){pow_CGM_95[2]=0};
pow_THO_95 = t(abs(beta_hat - true_beta - alt) > SE_THO*qnorm(0.975)); if( is.na(pow_THO_95)[1] ){pow_THO_95[1]=0}; if( is.na(pow_THO_95)[2] ){pow_THO_95[2]=0};
pow_DDG_95 = t(abs(beta_hat - true_beta - alt) > SE_DDG*qnorm(0.975)); if( is.na(pow_DDG_95)[1] ){pow_DDG_95[1]=0}; if( is.na(pow_DDG_95)[2] ){pow_DDG_95[2]=0};
pow_CHS_95 = t(abs(beta_hat - true_beta - alt) > SE_CHS*qnorm(0.975)); if( is.na(pow_CHS_95)[1] ){pow_CHS_95[1]=0}; if( is.na(pow_CHS_95)[2] ){pow_CHS_95[2]=0};

powers_IID_95 = cbind(powers_IID_95,pow_IID_95[1,2])
powers_CRi_95 = cbind(powers_CRi_95,pow_CRi_95[1,2])
powers_CRt_95 = cbind(powers_CRt_95,pow_CRt_95[1,2])
powers_CGM_95 = cbind(powers_CGM_95,pow_CGM_95[1,2])
powers_MNW_95 = cbind(powers_MNW_95,pow_MNW_95[1,2])
powers_THO_95 = cbind(powers_THO_95,pow_THO_95[1,2])
powers_DDG_95 = cbind(powers_DDG_95,pow_DDG_95[1,2])
powers_CHS_95 = cbind(powers_CHS_95,pow_CHS_95[1,2])
}

list_IID_95 = rbind(list_IID_95, powers_IID_95)
list_CRi_95 = rbind(list_CRi_95, powers_CRi_95)
list_CRt_95 = rbind(list_CRt_95, powers_CRt_95)
list_CGM_95 = rbind(list_CGM_95, powers_CGM_95)
list_MNW_95 = rbind(list_MNW_95, powers_MNW_95)
list_THO_95 = rbind(list_THO_95, powers_THO_95)
list_DDG_95 = rbind(list_DDG_95, powers_DDG_95)
list_CHS_95 = rbind(list_CHS_95, powers_CHS_95)

par(mar = c(5, 5, 2, 2))
plot(true_beta[2]+list_alt,colMeans(list_IID_95),type="l",lty=3,lwd=1,col="darkgray",ylim=c(0,1),xlab=expression(beta[1]),ylab="Power");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CRi_95),type="l",lty=2,lwd=1,col="darkgray",ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CRt_95),type="l",lty=1,lwd=1,col="darkgray",ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CGM_95),type="l",lty=3,lwd=1,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_MNW_95),type="l",lty=2,lwd=1,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_MEN_95),type="l",lty=1,lwd=1,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_THO_95),type="l",lty=2,lwd=2,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CHS_95),type="l",lty=1,lwd=2,ylim=c(0,1),xlab="",ylab="");
abline(v=true_beta[2],col="gray");abline(h=0.05,col="gray");
legend(min(true_beta[2]+list_alt),0.45,c("EHW","CRi","CRt","CGM","MNW","M","T","CHS"),lty=rev(c(1,2,1,2,3,1,2,3)),lwd=c(1,1,1,1,1,1,2,2),col=c("darkgray","darkgray","darkgray","black","black","black","black","black"))

options(digits=3)
print(c("Iteration",iter))
Sys.sleep(0.00000000000000001); flush.console()
}
########################################################################
# END OF MONTE CARLO ###################################################
########################################################################


par(mar = c(5, 5, 2, 2))
plot(true_beta[2]+list_alt,colMeans(list_IID_95),type="l",lty=3,lwd=1,col="darkgray",ylim=c(0,1),xlab=expression(beta[1]),ylab="Power");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CRi_95),type="l",lty=2,lwd=1,col="darkgray",ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CRt_95),type="l",lty=1,lwd=1,col="darkgray",ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CGM_95),type="l",lty=3,lwd=1,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_MNW_95),type="l",lty=2,lwd=1,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_MEN_95),type="l",lty=1,lwd=1,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_THO_95),type="l",lty=2,lwd=2,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CHS_95),type="l",lty=1,lwd=2,ylim=c(0,1),xlab="",ylab="");
abline(v=true_beta[2],col="gray");abline(h=0.05,col="gray");
legend(min(true_beta[2]+list_alt),0.45,c("EHW","CRi","CRt","CGM","MNW","M","T","CHS"),lty=rev(c(1,2,1,2,3,1,2,3)),lwd=c(1,1,1,1,1,1,2,2),col=c("darkgray","darkgray","darkgray","black","black","black","black","black"))

if(DGP_AR){
 filename = paste("results_AR_f_",b[1],"_",b[2],"_",b[3],"_rho",rho,"_N",N,"_T",T,".jpg",sep="")
}else{
 filename = paste("results_MA_f_",b[1],"_",b[2],"_",b[3],"_q",q,"_N",N,"_T",T,".jpg",sep="")
}
jpeg(filename, quality = 200)
par(mar = c(5, 5, 2, 2))
plot(true_beta[2]+list_alt,colMeans(list_IID_95),type="l",lty=3,lwd=1,col="darkgray",ylim=c(0,1),xlab=expression(beta[1]),ylab="Power");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CRi_95),type="l",lty=2,lwd=1,col="darkgray",ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CRt_95),type="l",lty=1,lwd=1,col="darkgray",ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CGM_95),type="l",lty=3,lwd=1,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_MNW_95),type="l",lty=2,lwd=1,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_MEN_95),type="l",lty=1,lwd=1,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_THO_95),type="l",lty=2,lwd=2,ylim=c(0,1),xlab="",ylab="");par(new=TRUE);
plot(true_beta[2]+list_alt,colMeans(list_CHS_95),type="l",lty=1,lwd=2,ylim=c(0,1),xlab="",ylab="");
abline(v=true_beta[2],col="gray");abline(h=0.05,col="gray");
legend(min(true_beta[2]+list_alt),0.45,c("EHW","CRi","CRt","CGM","MNW","M","T","CHS"),lty=rev(c(1,2,1,2,3,1,2,3)),lwd=c(1,1,1,1,1,1,2,2),col=c("darkgray","darkgray","darkgray","black","black","black","black","black"))
dev.off()
